1 데이터와 예측모형12
앞서 작업한 데이터와 선형회귀모형과 ARIMA 모형을 저장한 data/corona_fcst_list.rds 리스트 객체를 단변량 시계열 분석 작업을 위해서 풀어둔다.
library(tidyverse)
library(tidymodels)
library(timetk)
library(modeltime)
# 데이터 + 모형 ----
corona_fcst_list <- read_rds("data/corona_fcst_list.rds")
# 데이터 ----
full_tbl <- corona_fcst_list$data
# 모형 ----
wkfl_fit_lm <- corona_fcst_list$model$wkfl_fit_lm
wkfl_fit_arima <- corona_fcst_list$model$wkfl_fit_arima1.1 훈련/시험 데이터 분할
예측할 데이터와 모형개발에 활용할 데이터로 나눈 후에 모형개발에 활용할 데이터를 훈련/시험 데이터로 나눈다.
## 예측 데이터와 모형개발 데이터로 분리
forecast_tbl <- full_tbl %>%
filter(is.na(확진자))
history_tbl <- full_tbl %>%
filter(!is.na(확진자))
## 훈련/시험 데이터 분할
splits <- history_tbl %>%
time_series_split(date_var = 날짜,
assess = 30,
cumulative = TRUE)
splits<Analysis/Assess/Total>
<309/30/339>훈련/시험 데이터로 잘 나눠졌는지 시각적으로 확인한다.
splits %>%
tk_time_series_cv_plan() %>%
plot_time_series_cv_plan(.date_var = 날짜,
.value = 확진자)2 피처 공학
tidymodels 생태계의 recipe 팩키지 recipe() 함수를 사용해서 피쳐 공학(feature engineering) 작업을 수행한다. 추후, 다양한 피처를 작업할 예정이라 대략적인 틀만 잡아둔다.
recipe_spec <- recipes::recipe(확진자 ~ 날짜, data = training(splits))
recipe_spec %>% prep() %>% juice()# A tibble: 309 x 2
날짜 확진자
<date> <int>
1 2020-01-23 0
2 2020-01-24 1
3 2020-01-25 0
4 2020-01-26 1
5 2020-01-27 1
6 2020-01-28 0
7 2020-01-29 0
8 2020-01-30 0
9 2020-01-31 7
10 2020-02-01 1
# ... with 299 more rows3 단변량 예측 알고리즘
단변량 예측 알고리즘으로 다음 세가지 대표적인 알고리즘을 적용시킨다.
- ARIMA
- ETS
- TBATS
- Prophet
3.1 ARIMA
대표적인 시계열 데이터 예측모형 ARIMA로 적합시킨다.
mdl_arima_spec <- arima_reg() %>%
set_engine("auto_arima")
wkfl_fit_arima <- workflows::workflow() %>%
add_recipe(recipe_spec) %>%
add_model(mdl_arima_spec) %>%
fit(training(splits))
wkfl_fit_arima== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: arima_reg()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
Series: outcome
ARIMA(0,1,1)
Coefficients:
ma1
-0.2198
s.e. 0.0519
sigma^2 estimated as 2729: log likelihood=-1654.94
AIC=3313.87 AICc=3313.91 BIC=3321.333.2 ETS
지수평활(exponential smoothing) 중 ETS(Error-Trend-Season)를 예측모형으로 적합시킨다.
mdl_ets_spec <- exp_smoothing(
error = "auto",
trend = "auto",
season = "auto",
damping = "auto"
) %>%
set_engine("ets")
wkfl_fit_ets <- workflows::workflow() %>%
add_recipe(recipe_spec) %>%
add_model(mdl_ets_spec) %>%
fit(training(splits))
wkfl_fit_ets== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: exp_smoothing()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
ETS(A,N,A)
Call:
forecast::ets(y = outcome, model = model_ets, damped = damping_ets)
Smoothing parameters:
alpha = 0.7826
gamma = 1e-04
Initial states:
l = 2.4739
s = 1.6161 1.1827 -17.5283 -11.5425 14.7329 4.8602
6.6788
sigma: 51.3469
AIC AICc BIC
4216.526 4217.265 4253.860 3.3 TBATS
TBATS를 예측모형으로 적합시킨다.
mdl_tbats_spec <- seasonal_reg(
mode = "regression",
seasonal_period_1 = "auto"
) %>%
set_engine("tbats")
wkfl_fit_tbats <- workflows::workflow() %>%
add_recipe(recipe_spec) %>%
add_model(mdl_tbats_spec) %>%
fit(training(splits))
wkfl_fit_tbats== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: seasonal_reg()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
TBATS(1, {0,0}, -, {<7,2>})
Call: forecast::tbats(y = outcome, use.parallel = use.parallel)
Parameters
Alpha: 0.7812103
Gamma-1 Values: -0.002287667
Gamma-2 Values: 0.00209979
Seed States:
[,1]
[1,] 21.234291
[2,] 11.179662
[3,] -6.490442
[4,] 2.899411
[5,] -3.565564
Sigma: 50.80503
AIC: 4215.1043.4 Prophet
페이스북에서 개발한 Prophe 예측모형으로 적합시킨다.
mdl_prophet_spec <- prophet_reg(
) %>%
set_engine("prophet")
wkfl_fit_prophet <- workflows::workflow() %>%
add_recipe(recipe_spec) %>%
add_model(mdl_prophet_spec) %>%
fit(training(splits))
wkfl_fit_prophet== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: prophet_reg()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
PROPHET Model
- growth: 'linear'
- n.changepoints: 25
- changepoint.range: 0.8
- yearly.seasonality: 'auto'
- weekly.seasonality: 'auto'
- daily.seasonality: 'auto'
- seasonality.mode: 'additive'
- changepoint.prior.scale: 0.05
- seasonality.prior.scale: 10
- holidays.prior.scale: 10
- logistic_cap: NULL
- logistic_floor: NULL
- extra_regressors: 04 모형 평가
workflow 객체를 modeltime 팩키지 modeltime_table() 함수에 넣어 객체로 만든 후에, 추가로 지수평활법 EST 예측모형을 추가한다. modeltime_accuracy() 함수를 사용해서 모형 성능을 파악한다.
model_tbl <- modeltime_table(
wkfl_fit_arima,
wkfl_fit_ets,
wkfl_fit_tbats,
wkfl_fit_prophet
) %>%
update_model_description(.model_id = 1, "ARIMA") %>%
update_model_description(.model_id = 2, "ETS") %>%
update_model_description(.model_id = 3, "TBATS") %>%
update_model_description(.model_id = 4, "Prophet")
model_tbl %>%
modeltime::modeltime_accuracy(testing(splits))# A tibble: 4 x 9
.model_id .model_desc .type mae mape mase smape rmse rsq
<int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 ARIMA Test 275. 30.3 2.78 36.7 341. NA
2 2 ETS Test 277. 30.5 2.80 37.0 343. 0.0720
3 3 TBATS Test 277. 30.5 2.81 37.1 343. 0.0650
4 4 Prophet Test 603. 73.0 6.10 116. 647. 0.461 4.1 시각화
시험데이터를 통해 모형으로 예측한 값을 시각화한다. ETS MAE 값이 ARIMA 예측모형과 별다른 차이가 나고 있지 않지만, 신뢰구간이 생긴 것은 그나마 다행이다.
calibration_tbl <- model_tbl %>%
modeltime_calibrate(
new_data = testing(splits)
)
calibration_tbl %>%
modeltime_forecast(
new_data = testing(splits),
actual_data = history_tbl,
keep_data = TRUE,
conf_interval = 0.10
) %>%
plot_modeltime_forecast(
.legend_max_width = 60,
.legend_show = TRUE,
.conf_interval_show = TRUE,
.conf_interval_alpha = 0.20,
.conf_interval_fill = "lightblue",
.title = "코로나19 확진자 1개월 예측"
)5 확진자 예측
앞선 모형은 history_tbl 을 훈련/시험 데이터에 대해 적합을 시킨 것이라 … 이제 시간을 확대하여 modeltime_refit() 함수를 사용해서 모형을 전체 데이터에 대해 다시 적합시킨다.
refit_tbl <- calibration_tbl %>%
modeltime_refit(data = history_tbl)마지막으로 앞서 구축된 모형을 바탕으로 현재까지 입수된 데이터를 바탕으로 한달 후를 예측한다.
refit_tbl %>%
modeltime_forecast(
new_data = forecast_tbl,
actual_data = history_tbl,
conf_interval = 0.3
) %>%
plot_modeltime_forecast(
.legend_max_width = 25,
.conf_interval_fill = "lightblue",
.conf_interval_alpha = 0.7,
.interactive = TRUE
)6 모형과 데이터 저장
마지막 단계로 데이터와 모형을 저장하여 다음 단계로 나가기 위한 작업을 수행한다.
corona_fcst_list <- list(
# 데이터 ---
data = full_tbl,
# 예측모형 ----
model = list(
wkfl_fit_lm = wkfl_fit_lm,
wkfl_fit_arima = wkfl_fit_arima,
wkfl_fit_ets = wkfl_fit_ets,
wkfl_fit_tbats = wkfl_fit_tbats,
wkfl_fit_prophet = wkfl_fit_prophet
)
)
# 모형 + 데이터 저장 ----
corona_fcst_list %>%
write_rds("data/corona_fcst_list.rds")데이터 과학자 이광춘 저작
kwangchun.lee.7@gmail.com